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Stirring  and  mixing  of  biologically  reactive  tracers 

K.  J.  Richards,  S.  J.  Brentnall,  R McLeod,  and  A.  P.  Martin 

Southampton  Oceanography  Centre,  Southampton,  UK 


Abstract.  We  examine  the  effect  of  stirring  and  mixing  on  the  marine 
planktonic  ecosystem.  We  present  a number  of  models  investigating  the 
impact  of  patchiness  in  the  population,  diffusion,  and  fluid  stirring  on  the 
spatial  structure  of  the  population  and  the  overall  production  rates  of  the 
system.  It  is  found  that  the  heterogeneity  of  the  system  significantly  affects 
the  way  the  system  behaves.  The  response  of  the  ecosystem  to  the  action  of 
the  fluid  flow  is  a function  of  not  only  the  topology  of  the  flow  but  also  the 
dynamics  of  the  ecological  model  itself. 


1.  Introduction 

The  marine  ecosystem  plays  an  important  role  in  the 
cycling  of  carbon.  In  particular  it  acts  as  a ‘carbon 
pump,’  utilising  atmospheric  CO2  through  photosyn- 
thesis, a proportion  of  the  carbon  being  exported  to 
the  deep  ocean,  and  thus  moderating  levels  of  CO  2 in 
the  atmosphere.  Ecological  models  are  now  incorpo- 
rated into  models  of  climate  change  to  predict  global 
warming  scenarios  (see,  e.g.,  it  Cox  et  al.,  2000).  A 
physical/ecological  model  used  in  climate  studies  needs 
to  incorporate  both  the  appropriate  ecological  dynam- 
ics and  the  right  physics  impacting  on  the  biological 
production. 

Primary  production  in  the  ocean  is  dominated  by 
phytoplankton.  A feature  of  the  distribution  of  phyto- 
plankton distributions  is  that  it  is  very  heterogeneous 
or  ‘patchy’  ( Bainbridge , 1957).  Structure  is  observed 
on  scales  ranging  from  metres  to  the  basin  scale  ( Mann 
and  Lazier , 1996)  with  this  structure  often  observed  to 
be  associated  with  physical  features  such  as  eddies  and 
fronts  (e.g.,  (The  Ring  Group,  1981;  Falkowski  et  al., 
1991;  Strass,  1992).  A number  of  important  questions 
arise: 


• Does  the  patchiness  affect  the  dynamics  of  the 
ecological  system? 

• Does  the  heterogeneity  of  the  system  impact  on 
the  overall  production  rates  and/or  community 
structure? 

• Does  the  structure  affect  the  way  the  biological 
system  responds  to  changes  to  the  environment? 


damental  issues  related  to  the  problem  of  marine  bio- 
logical patchiness.  The  models  should  be  thought  of  as 
caricatures  of  the  more  complex  system,  but  we  hope 
they  will  provide  the  building  blocks  for  -understanding 
of  that  system.  The  work  is  very  much  ‘in  progress’ 
and  so  this  not  a presentation  of  a ‘complete  theory’ 
but  rather  some  initial  attempts  and  their  results. 


2.  Biological  reactions  and  fluid  flow 

Fluid  flow  can  impact  on  the  marine  ecological  sys- 
tem in  two  distinct  ways:  (a)  through  the  vertical  move- 
ment of  nutrients  and  biological  species  in  or  out  of 
the  euphotic  zone,  and  (b)  the  stirring  and  mixing  of 
components  of  the  system.  A good  example  of  when 
both  are  important  is  baroclinic  instability  ( Spall  and 
Richards,  2000).  Associated  with  the  instability  process 
is  a permanent  vertical  displacement  carrying  nutrient 
rich  waters  into  the  sunlit  surface  layers.  These  nutri- 
ents are  then  stirred  and  mixed  with  the  surrounding 
water  and  are  utilised  by  the  biota.  Here  we  will  con- 
centrate on  this  second  process,  the  stirring  and  mixing 
by  lateral  fluid  flow. 


2.1.  General  advection-diffusion-reaction 
equation 

We  shall  model  the  ecological  system  using  the  advec- 
tion-diffusion-reaction equation 

dE 

— -f  (u  - V)  E = f (E)  + V-  (DVE) 
where,  for  example,  the  state  vector  may  take  the  form 


E = 


N 

P 

Z 


Here  we  will  use  a number  of  relatively  simple  mod- 
els of  both  the  ecological  dynamics  and  the  fluid  flow 
in  order  to  try  and  understand  some  of  the  more  fun- 
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with  the  elements  representing  the  concentration  of  nu- 
trients (IV),  phytoplankton  (P)  and  zooplankton  ( Z ) 
expressed  in  terms  of  a currency,  such  as  nitrogen. 

The  elements  are  advected  by  the  flow  u and  dif- 
fused. The  reaction  between  elements  is  represented  by 
f(E).(How  well  the  field  equation  approach  models  the 
marine  ecological  system  is  a matter  of  debate.  There 
is  some  evidence  that  it  fails  if  behavioural  aspects  of 
the  zooplankton  are  important  ( Flierl  et  al,  1999). 

A measure  of  the  relative  ‘strength’  of  the  reaction 
terms  to  either  diffusion  or  advection  can  be  gauged  by 
considering  the  ratio  of  timescales,  i.e. 

K 

and 


respectively,  where  U and  L are  typical  velocity  and 
length  scales,  k a diffusion  coefficient,  and  fi  the  reaction 
rate.  The  nondimensional  parameters  0K  and  3U  are 
equivalent  to  diffusive  and  advective  inverse  Damkohler 
numbers.  Large  values  imply  that  the  reaction  domi- 
nates. 

We  present  three  models  using  (2.1).  The  first  looks 
at  the  case  of  no  flow  and  the  structure  that  evolves  sim- 
ply through  the  combination  of  reaction  and  diffusion. 
The  second  considers  simple  biology  in  a simple  flow,  a 
pure  strain,  and  examines  the  structure  of  filaments  of 
biological  material.  The  third  takes  a 2D  geostrophic 
turbulent  flow  and  investigates  the  impact  the  stirring 
has  on  biological  production.  The  underlying  message 
is  that  the  response  of  the  ecology  to  the  action  of  the 
fluid  flow  is  a function  of  not  only  the  topology  of  the 
flow  but  also  the  dynamics  of  the  ecological  model  itself. 


as  Holling  type  III).  The  zooplankton,  Z,  grow  at  a re- 
duced rate  7 < 1.  The  mortality  term  h(Z)  is  taken 
to  be  either  a linear  or  quadratic  function  of  Z,  corre- 
sponding to  the  zooplankton  dying  of  old  age  or  being 
consumed  by  a higher  predator,  respectively. 
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Figure  1.  Model  1 with  linear  mortality.  Upper  panel, 
plot  of  phytoplankton  concentration  as  a function  of  distance 
(given  in  tens  of  metres)  and  time.  Lower  panel,  trajectories 
of  the  mean  values  of  P and  Z for  the  randomly  perturbed 
(crosses)  and  unperturbed  (dots)  cases.  Also  shown  are  the 
null-clines  of  the  reaction  equations  (i.e.,  dP/dt  — 0 and 
dZ/dt  = 0.  The  null-clines  intersect  giving  a single  stable 
equilibrium  point 


2.2.  Model  1:  Intrinsic  structure 

Consider  the  case  of  no  flow  (u  = 0).  We  shall  con- 
sider a two  compartment  model: 


where  P € [0, 1],  and  the  reaction  terms  are 


The  model  for  the  reaction  terms  is  that  of  Truscott 
and  Brindley  (1994).  It  is  chosen  because  of  its  sim- 
plicity and  because  the  mathematical  structure  of  the 
ecological  dynamics  is  known.  In  particular  the  system 
becomes  excitable  (i.e.,  small  but  finite  perturbations 
from  an  equilibrium  point  can  lead  to  large  excursions 
in  model  phase  space)  if  both  7 and  v are  small  enough. 

The  diffusion  tensor  is  taken  to  be 


f = 


rP(l-P)-&* 

7 (zrrp*)  ~uh{Z) 


The  reaction  terms  for  the  phytoplankton,  P,  include 
limited  growth  (known  as  logistic  growth)  and  grazing 
by  zooplankton  (the  functional  form  used  here  is  known 


so  that  zooplankton  may  diffuse  at  a different  rate  to 
phytoplankton  (allowing  cross  diffusion,  i.e.,  zooplank- 
ton diffusing  up  the  phytoplankton  gradient,  has  also 
been  investigated). 
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Here  we  present  two  examples  to  illustrate  the  dif- 
ferent behaviour  of  the  system  for  linear  (Figure  1)  and 
quadratic  (Figure  2)  mortality  (further  details  of  the 
behaviour  of  the  system  are  given  in  Brentnall  et  d. 
(2001).  We  take  a one-dimensional  case,  i.e.,  variables 
only  vary  in  one  horizontal  direction.  The  system  is  ini- 
tialised on  a grid  with  random  values  of  P and  Z about 
some  mean  value.  The  size  of  the  rectangular  box  in 
the  lower  panel  of  each  figure  gives  an  indication  of  the 
spread  of  values.  In  both  the  linear  and  quadratic  mor- 
tality cases  a structure  appears.  The  difference  is  that 
this  structure  persists  for  a longer  time  for  the  quadratic 
mortality  case.  The  length  scale  of  the  structure  scales 
with  y/n/n  (i.e.,  the  diffusive  inverse  Damkohler  num- 
ber, j3K  is  the  relevant  parameter).  With  a growth  rate 
of  1 day-1  and  diffusion  coefficient  of  10~2  m2s_1  the 
emergent  scale  is  of  order  ~ 500  m. 
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Figure  2.  Model  1 with  quadratic  mortality.  As  for  Figure 
1,  but  in  this  case  the  null-clines  intersect  at  three  points, 
with  two  stable  equilibria  at  P — 0.05  and  0.7  and  an  un- 
stable equilibrium  point  at  P = 0.27 

The  trajectory  in  P/Z  phase  space  changes  signif- 
icantly when  horizontal  structure  is  introduced.  The 
trajectories  for  the  mean  P and  Z with  and  without 
the  small  random  perturbation  are  shown  in  the  lower 
panel  of  Figures  1 and  2.  With  linear  mortality,  al- 
though the  two  trajectories  converge  to  the  same  stable 
equilibrium,  the  excursion  in  P is  much  greater  when 


the  perturbation  is  introduced.  With  quadratic  mortal- 
ity, the  two  trajectories  converge  to  two  different  sta- 
ble equilibria,  with  the  perturbed  case  converging  very 
slowly.  Clearly  the  heterogeneity  of  the  biology  (with 
diffusion  acting)  is  playing  a significant  role  in  the  time 
evolution  of  the  system.  It  is  also  clear  that  we  need  to 
consider  what  happens  when  there  is  a flow,  which  we 
do  in  the  next  section 

2.3.  Model  2:  filaments 

We  now  consider  a pure  strain  flow 

u = ( -Xx , A y) 

Such  a flow  will  tend  to  pull  a tracer  into  a thin  filament. 
For  an  inert  tracer  the  filament  width  attains  a value 


when  the  straining  of  the  flow  is  balanced  by  diffusion. 
(34  (@)  has  shown  that  filaments  of  an  exponentially 
growing  tracer  have  the  same  width  as  an  inert  tracer) 

The  biology  we  will  use  is  a simplified  version  of  that 
used  above,  namely 

E=[P] 

with 

f=  [ /iP(l-P)  ] 

A feature  of  this  system  is  that  the  combination  of 
logistic  growth  and  diffusion  allows  the  existence  of  re- 
active travelling  waves  (or  Fisher  waves)  which  have 
a minimum  speed  2v/7cp  (c.f.,  Murray,  1993).  These 
waves  are  important  in  transmitting  ‘information’  and 
means  that  the  biological  population  can  advance  faster 
than  through  pure  diffusion.  The  waves  are  important 
in  setting  the  width  of  filaments. 

We  may  expect  the  travelling  wave  to  be  brought  to 
rest  when  its  speed  matches  that  of  the  converging  flow, 
such  that  the  -width  of  the  filament 


Indeed  this  is  found  to  be  the  case  for  sufficiently 
large  values  of  (3U  = fi/X,  an  inverse  advective  Damkohler 
number  (see  the  lower  panel  of  Figure  3).  The  shape  of 
the  distribution  of  concentration  changes  with  increas- 
ing j3u  (upper  panel  Figure  3)  from  a Gaussian  to  a 
square  wave  as  shown  by  the  kurtosis  of  the  distribu- 
tion (lower  panel  Figure  3)  ( McLeod  et  al.,  2001). 
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Figure  3.  Model  2,  filament  structure.  Upper  panel:  the 
distribution  of  phytoplankton  across  the  filament  for  three 
values  of  /?„:  dash-dot  /3U=2,  dash  0U— 10  and  solid  /3U=50. 
The  horizontal  distance  has  been  scaled  with  y/ Kfl/X . Lower 
panel:  the  filament  width  (scaled  with  y/n/X)  and  kurtosis 
as  a function  of  0U. 


potential  vortlctty 


Figure  4.  Model  3,  fluid  stirring.  The  potential  vorticity 
field. 


The  biological  model  is  slightly  more  complicated 
than  before. 


E = 


N 

P 

Z 

D 


with 


It  is  interesting  to  note  that 


Wr 

Wi 


V~0 


i.e.,  the  filament  width  for  a reactive  tracer  will  be 
greater  than  that  of  an  inert  tracer  for  (3U  > 1.  Using  an 
observed  filament  width,  and  strain  rate,  to  estimate  k 
will  lead  to  an  overestimate  by  a factor  0 if  the  reaction 
is  not  taken  into  account. 


2.4.  Model  3:  fluid  stirring 

Lastly  we  consider  the  case  of  the  ecosystem  be- 
ing stirred  by  a turbulent  flow  field.  We  use  the  2D 
geostrophic  turbulence  model  of  Babiano  et  al.  (1987) 
in  a 512  km2  periodic  domain.  The  flow  is  forced  to 
give  a statistically  steady  state.  An  example  of  the  po- 
tential vorticity  is  given  in  Figure  4.  A characteristic 
of  2D  turbulence  is  that  the  flow  is  inhabited  by  strong 
coherent  vortical  structures  surrounded  by  a straining 
flow  in  the  intervening  regions.  The  coherent  structures 
act  as  transport  barriers  inhibiting  the  mixing  of  tracers 
between  the  eddies  and  surrounding  fluid  ( Provenzale , 
1999). 


-/*  (tACTn)  P + + 72z 

»(1^n)p-(1^)z-vpP 

Z-l2Z-vzZ* 

(1  - 7i)  (0x)  Z + ixpP  + ~ VdD 


where  D represents  detritus  and  there  is  an  additional 
sinking  term  —wsD/h  in  the  equation  for  detritus.  The 
biological  model  is  that  of  Oschlies  and  Garcon  (1999) 
and  has  been  chosen  because  it  has  been  shown  to  ca- 
pable of  reproducing  the  seasonal  cycle  in  both  olig- 
otrophic  (nutrient  depleted)  and  nonoligotrophic  re- 
gions with  the  same  set  of  parameter  values.  Details 
of  the  biological  and  flow  model,  including  the  parame- 
ter values,  are  given  in  Martin  et  al.  (2001). 

The  biological  system  is  forced  by  a source  of  nutrient 


rs  - s{N0  - N) 


where  s is  a function  of  position  or  some  property  of  the 
flow,  and  can  be  thought  of  as  modelling  the  input  of  nu- 
trients to  the  euphotic  zone  through  upwelling.  Martin 
et  al.  (2001)  investigate  how  the  productivity  depends 
upon  the  distribution  of  the  upwelling  by  comparing 
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spatially  stationary  sources  with  sources  correlated  or 
decorrelated  with  flow  features.  They  conclude  that  the 
increased  productivity  brought  about  by  upwelling  is  a 
function  of  the  rate  of  upwelling,  the  distribution  of  the 
sources,  and  the  mixing  efficiency  of  the  flow. 

We  consider  two  cases,  the  first  with  the  source  con- 
fined to  a single  patch  (Figure  5)  and  the  second  with 
the  source  distributed  over  64  smaller  patches  but  such 
that  the  total  area  is  the  same  as  in  the  first  case  (Fig- 
ure 6).  The  figures  show  the  distribution  of  primary 
production  at  a time  when  the  biology  is  in  a statistical 
equilibrium. 


forcing 


primary  production  (mMol  N rrT3  d ’) 


Figure  5.  Model  3,  fluid  stirring.  The  forcing  of  nutrient 
and  primary  production  resulting  from  stirring  by  a 2D  tur- 
bulent field.  The  production  ranges  from  0.1  (black)  to  0.45 
(white)  nMol  N m~3  d-1. 

Referring  to  the  single  source  case  (Figure  5),  the 
nutrient  forcing  has  produced  elevated  rates  of  produc- 
tion above  the  background  level  of  0.1  mMol  N m~3 
d-1.  Compared  to  the  case  with  the  same  forcing  but 
no  fluid  stirring  (u  = 0)  the  areal  average  of  produc- 


tion has  increased  by  36%.  The  fluid  motion  has  the 
effect  of  increasing  the  surface  area  between  the  nutri- 
ent enriched  and  surrounding  waters.  A striking  feature 
of  the  distribution  of  production  is  that  it  is  limited  to 
the  vicinity  of  the  forcing  (an  inert  tracer  will  be  spread 
quickly  across  the  entire  domain).  The  reaction  limits 
the  extent  of  the  region  of  high  nutrients  and  hence  the 
area  of  contact  with  the  surrounding  waters.  Reduc- 
ing the  reaction  rate  (reducing  (3U)  will  allow  greater 
stirring  to  take  place  before  the  nutrients  are  consumed 
and  hence  increase  the  area  of  contact,  counteracting 
the  decrease  in  overall  production  because  of  the  de- 
creased production  rate.  Experiments  varying  @u  show 
that  the  two  effects  may  balance  each  other. 


forcing 
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km 


Figure  6.  As  for  Figure  5 but  for  a distributed  source. 


As  expected,  distributing  the  source  increases  the 
mixing  and  hence  production  rate  (Figure  6).  In  this 
case  the  overall  production  rate  has  been  increased  by 
137%  above  the  case  with  no  mixing. 
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3.  Concluding  remarks 

We  have  reported  on  a number  of  investigations  into 
the  impact  of  heterogeneity  of  the  marine  ecological  sys- 
tem on  the  behaviour  of  that  system  as  a whole.  The 
results  demonstrate  that  the  heterogeneity  cannot  be  ig- 
nored. Applying  the  results  from  studies  such  as  those 
presented  here  is  difficult.  We  freely  admit  that  the 
results  of  our  studies  to  date  axe  somewhat  of  an  anec- 
dotal character.  There  are  caveats  on  the  methodology 
such  as  the  exclusion  of  the  diurnal  and  in  particular 
the  seasonal  cycle,  both  of  which  impose  a strong  pac- 
ing of  the  system,  and  which  need  to  be  included  in 
future  studies. 

It  is  disconcerting  that  the  results  are  very  depen- 
dent on  the  form  of  the  ecological  model  (as  an  addi- 
tional example,  the  enhancement  of  productivity  found 
in  model  3 was  significantly  reduced  when  a simpler 
[N,P]  was  used).  Perhaps  not  a surprising  remark  but 
one  which  is  often  ignored.  This  issue  needs  to  be  ex- 
plored further,  but  it  does  make  the  task  of  developing 
‘robust’  ecological  models  that  much  more  difficult. 

What  is  clear  is  that  viewing  marine  ecology  as  a 
simple  ID  system  (the  space  dimension  being  vertical) 
may  well  be  erroneous  and  produce  misleading  results. 
Ecological  models  are  often  ‘fitted’  or  rejected  on  the 
basis  of  the  performance  of  a ID  version  of  the  model 
compared  to  observations.  At  best  the  ‘fitted’  param- 
eters may  be  dependent  on  the  physical  environment, 
implicitly  including  the  effects  of  unresolved  physical 
and  ecological  processes.  At  worst  the  functional  form 
of  the  model  chosen  through  this  comparison  may  be 
wrong.  In  either  case,  in  using  an  ecological  model  to 
predict  changes  to  the  ecological  system  as  the  physical 
environment  changes,  we  need  to  ensure  that  the  model 
captures  the  correct  impact  of  those  changes. 
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